function ground_reaction_force(nlp, bounds)
domain = nlp.Plant;
f_right = SymVariable('fright',[3,1]);
f_left = SymVariable('fleft',[3,1]);
GRF = [f_right(3);f_left(3)];
GRF_func = SymFunction(['GroundReactionForce_',domain.Name],GRF,{f_right,f_left});
if strcmp(domain.Name,'FrontStance')
    GRF_cstr = NlpFunction('Name',GRF_func.Name,...
        'Dimension',2,...
        'lb', bounds.constrBounds.GRF.lb,...
        'ub', bounds.constrBounds.GRF.ub,...
        'Type', 'Linear',...
        'SymFun', GRF_func, ...
        'DepVariables',[nlp.OptVarTable.fRightFront(1);nlp.OptVarTable.fLeftFront(1)]);
    addConstraint(nlp,GRF_cstr.Name,'first',GRF_cstr);
else
    GRF_cstr = NlpFunction('Name',GRF_func.Name,...
        'Dimension',2,...
        'lb', bounds.constrBounds.GRF.lb,...
        'ub', bounds.constrBounds.GRF.ub,...
        'Type', 'Linear',...
        'SymFun', GRF_func, ...
        'DepVariables',[nlp.OptVarTable.fRightBack(1);nlp.OptVarTable.fLeftBack(1)]);
    addConstraint(nlp,GRF_cstr.Name,'first',GRF_cstr);
end
end

